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We study the non-linear evolution of baryon acoustic oscillations in the matter power spectrum 
and correlation function from the improved perturbation theory (PT). Based on the framework of 
renormalized PT, which provides a non-perturbative way to treat the gravitational clustering of 
large-scale structure, we apply the closure approximation that truncates the infinite series of loop 
contributions at one-loop order, and obtain a closed set of integral equations for power spectrum 
and non-linear propagator. The resultant integral expressions are basically equivalent to those 
previously derived in the form of evolution equations, and they keep important non-perturbative 
properties which can dramatically improve the prediction of non-linear power spectrum. Employing 
the Born approximation, we then derive the analytic expressions for non-linear power spectrum and 
the predictions are made for non-linear evolution of baryon acoustic oscillations in power spectrum 
and correlation function. We find that the improved PT possesses a better convergence property 
compared with standard PT calculation. A detailed comparison between improved PT results 
and N-body simulations shows that a percent-level agreement is achieved in a certain range in 
power spectrum and in a rather wider range in correlation function. Combining a model of non- 
linear redshift-space distortion, we also evaluate the power spectrum and correlation function in 
redshift space. In contrast to the results in real space, the agreement between N-body simulations 
and improved PT predictions tends to be worse, and a more elaborate modeling for redshift-space 
distortion needs to be developed. Nevertheless, with currently existing model, we find that the 
prediction of correlation function has a sufficient accuracy compared with the cosmic- variance errors 
for future galaxy surveys with volume of a few h~^Gpc^ at z > 0.5. 

PACS numbers: 98.80.-k 



I. INTRODUCTION 



In the last decade, systematic measurements of the cosmic microwave background anisotropies as well as large-scale 
structure of the Universe have led to the estabhshment of the "standard cosmological model" (e.g., [l|,[3,[l,l3,IM)- The 
Universe is close to a flat geometry, and is filled with the hypothetical cold dark matter (CDM) particles, together 
with a small fraction of baryons, which serve as the seeds of structure formation of the Universe. The most striking 
feature in the standard cosmological model is that the energy contents of the Universe is dominated by the mysterious 
energy component called dark energy, which is supposed to drive the late-time cosmic acceleration discovered by the 
observation of distant supernovae (e.g., [1, 0|)- 

Currently, our understanding of the nature of dark energy is still lacking. Although the observation is roughly 
consistent with cosmological constant and with no evidence for time dependence of dark energy, long-distance modifi- 
cations of general relativity have been proposed alternative to the dark energy and these reconcile with the observation 
of late-time acceleration (see [8l M, [To, 11] for reviews) . While a fully consistent model of modified gravity has not yet 
been constructed ( see [ij, for popular models), a possibility of break-down of general relativity still remains 

and should be tested. 

To understand deeply the nature of dark energy or origin of cosmic acceleration, a further observational study 
is definitely important. There are two comprehensive ways to distinguish between many models of dark energy and 
discriminate the dark energy from modified gravity. One is to precisely measure the expansion history of the Universe, 
and the other is to observe the growth of structure. 

Among various observational techniques, baryon acoustic oscillations (BAOs) imprinted on the matter power spec- 
trum or two-point correlation function can be used as a standard ruler to measure the cosmic expansion history (e.g.. 
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[Tg} , see also [13, [iM [3 HO, HH for recent BAO measurements) . The characteristic scale of BAOs, which is de- 
termined by the sound horizon scale of primeval baryon-photon fluid at the last scattering surface [H, [Hi , is thought 
to be a robust measure and it lies on the linear to quasi-linear regimes of the gravitational clustering of large-scale 
structure [13, [IH]- With a percent-level determination of the characteristic scale of BAOs, the expansion history can 
be tightly constrained, and the equation-of-state parameter of the dark energy, Wdc, defined by the ratio of pressure 
to energy density of dark energy, would be precisely determined within the precision of a few % level [26l . l27l | . This is 
the basic reason why most of the planned and ongoing galaxy redshift surveys aim at precisely measuring the BAOs 



While the robustness of the BAOs as a standard ruler has been repeatedly stated and emphasized in the literature, 
in order to pursue an order-of-magnitude improvement, a precise theoretical modeling of BAOs definitely plays an 
essential role for precision measurement of BAO scale, and it needs to be investigated taking account of the various 
systematic effects. Among these, the non-linear clustering and redshift-space distortion effects as well as the galaxy 
biasing cannot be neglected, and affect the characteristic scale, although their effects are basically moderate at the 
relevant wavenumber, k < 0.3/iMpc~^. 

Recently, several analytic approaches to deal with the non-linear clustering have been developed, complementary to 
the N-body simulations |32, JMi, iiS,, j36j, J38,, ji9 , 4^, |4l|, (4^ . In contrast to the standard analytical calculation with 
perturbation theory (PT), these have been formulated in a non-perturbative way with techniques resumming a class 
of infinite series of higher-order corrections in perturbative calculation. Thanks to its non-perturbativc formulation, 
the applicable range of the prediction is expected to be greatly improved, and the non-linear evolution of baryon 
acoustic oscillations would be accurately described with a percent-level precision. 

The purpose of this paper is to investigate the viability of this analytic approach, focusing on a specific improved 
treatment. In the previous paper [43j . we have applied a non-linear statistical method, which is widely accepted 
in the statistical theory of turbulence [3], to the cosmological perturbation theory of large-scale structure. We 
have derived the non-perturbative expressions for the power spectrum, coupled with non-linear propagator, which 
effectively contain the information on the infinite series of higher-order corrections in the standard PT expansion. 
Based on this formalism, the analytic treatment of the non-perturbative expression is developed employing the Born 
approximation, and the leading-order calculation of power spectrum is compared with N-body simulations in real 
space [1^, finding that a percent-level agreement is achieved in a mildly non-linear regime (see also [11]). Here, we 
extend the analysis to those including the next-to- leading order corrections of Born approximation. In addition to the 
power spectrum, we will consider the two-point correlation function, paying a special attention on the baryon acoustic 
peak, i.e., a Fourier counterpart of BAOs in power spectrum. Further, we also discuss the non- linear clustering in 
redshift space, and the predictions of improved PT are compared with N-body results, combining a non-linear model 
of redshift-space distortion. We examine how well the present non-linear model accurately describe the systematic 
effects on BAOs and/or baryon acoustic peak. 

This paper is organized as follows. In Sec. |lTl we briefly mention the basic equations for cosmological PT as our 
fundamental basis to deal with the non-linear gravitational clustering. We then discuss in some details in Sec. IIIII 
how to compute the non- linear power spectrum or two-point correlation functions. Starting from the discussions on 
standard treatment of perturbative calculation and its non-perturbative reformulation called renormalized PT, we 
introduce the closure approximation, which gives a consistent non-perturbative scheme to treat the infinite series of 
renormalized PT expansions, and obtain a closed set of non-perturbative expressions for power spectrum. Based on 
this, we present a perturbative treatment of the closed set of equations while keeping important non-perturbative 
properties. Section IIVI gives the main result of this paper, in which a detailed comparison between improved PT 
calculation and N-body simulation is made, especially focusing on the non-linear evolution of BAOs. We compute the 
power spectrum and two-point correlation function in both real and redshift spaces, and investigate the accuracy of 
both predictions by comparing improved PT with N-body results. Finally, section |V] is devoted to the discussion and 
conclusion. 



Throughout the paper, we consider the evolution of cold dark matter (CDM) plus baryon systems neglecting the tiny 
fraction of (massive) neutrinos. Owing to the single-stream approximation of the coUisionless Boltzmann equation, 
which is thought to be a quite accurate approximation on large scales, the evolution of the CDM plus baryon system 
can be treated as the irrotational and pressureless fluid system whose governing equations are continuity and Euler 
equations in addition to the Poisson equation (see Ref. |47| for review). In the Fourier representation, these equations 
are further reduced to a more compact form. Let us introduce the two-component vector fe.g..[33|): 



(e.g., 28, 29, 30, 31]). 



II. PRELIMINARIES 




(2.1) 
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where the subscript a = 1, 2 selects the density and the velocity components of CDM plus baryons, with S and 
9{x) = V • v{x)/{aH), where a and H are the scale factor of the Universe and the Hubble parameter, respectively. 
The function f{t) is given by f{t) = d\nD{t) / dhia and the quantity D{t) being the linear growth factor. Then, in 
terms of the new time variable 77 = lnZ)(t), the evolution equation for the vector quantity <&Q(fc; t) becomes 



d 



/Ski (Pk2 
— -^--3— - fci - fca) 7ahc(fci, ^2) $b(fci; 77) $c(fe2; 77), 



(2.2) 



where Sd is the Dirac delta function. Here and in what follows, we use the summation convention that the repetition 
of the same subscripts indicates the sum over the whole vector components. The time-dependent matrix flabiv) is 
given by 







-1 



(2.3) 



The quantity flaiiv) is the density parameter of CDM plus baryons at a given time. Each component of the vertex 
function 'jabc becomes 



7afcc(fcl,fc2) = < 



' H^ + fel^} = («'^c) = (l,l,2) 



(2.4) 



(fcj ■fc2) |fcl +^2 
2|fei|2|fe2P 



(a,5,c) = (2,2,2) 
; otherwise 

Note that the formal solution of <i>Q can be obtained from Eq. (|2.2p and is expressed as (e.g., [H, H^l) 

' 3 ^ Snik - fei - fc2) jbcd{ki,k2)Mki;v')Mk2; v')- (2.5) 



*a(fc; v) = gabiv, Vo) Ub So{k) + / dT]'gab{r], v') 

J Vo 



(27r)3 



Here, the quantity Ua is the constant vector which specifies the initial condition (see next section), and the quantity 
gab denotes the linear propagator satisfying the following equation: 



d 

Sab-^ + ^abiv) 

orj 



9bc{v,il') = 0, 



(2.6) 



with the boundary condition gabiv^v) = ^ab- The quantity Sq is the random density field given at an early time r/o, 
which is assumed to obey the Gaussian statistic. The power spectrum of density field is defined as 



(<5o(fc)<5o(fe')) - (2vr)3 Soik + k') P„{k). 



(2.7) 



Eq. (|2.2p or (|2.5p is the fundamental building block of large-scale structure, and the three quantities ^abc, gab and 
PgUaUf, introduced here constitute the basic pieces of standard PT. The graphical representation of them is shown in 
Fig. [2 (see also Ref. [13]). 



III. IMPROVED PERTURBATION THEORY 



A. Standard PT vs. Renormalized PT 



In this paper, we are especially concerned with the non-linear evolution of the two-point statistics, defined as the 
ensemble average of ^a- 

($a(fe; 77)*b(fc'; ??')) = (27r)3(5D(fc + fc')i"afc(|fc|;r?,r;'); V > v' ■ (3-1) 

In the above, there are four types of power spectra, Pn, P12, P21 and P22, which respectively correspond to the auto- 
and cross-power spectra, Pss, —Pse/f, —Pes/f and Pee/P- Note that in general we have P12 ^ P21 unless 77 = 77'. 
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FIG. 1; Diagrammatic notion of the initial power spectrum (left), linear propagator (middle), and tree vertex (right). The 
linear propagator satisfies the equation (|2.6p with boundary condition Qabiji, ??) = Sab- The explicit expression of vertex function 
7a6c is given by Eq. (|2.4p . 



Consider how to compute the power spectrum based on the analytic treatment. In the standard treatment of the 
perturbation theory, we first assume that the field $a is a small perturbed quantity and it is expanded as 

^a{k; v) - (fe; v) + {k; v) + {k;v) + --- ■ (3.2) 

The explicit functional form of the quantity is systematically derived through the order-by-order treatment of 
Eq. (|2.2p . Substituting the above expansion into the definition (|3.ip and evaluating it perturbatively, the power 
spectrum Paf,(fc; 77, 77), shortly abbreviated as Pab{k;ri), is schematically expressed as 

Pabik; v) = e^^uaub Po(fc) + V) + P!l'°°''{k; v) + ■ ■ ■ ■ (3.3) 

where we chose Ua — (1, 1), which implies that the growing-mode solution is imposed at the initial condition^. The 
function Pa{k) is the linear power spectrum given at an early time, obtained from the first-order quantity C'l^'' 

(see Eq. (|2.7I for definition). The subsequent terms P^b°°^ and P^t°°^ represent the corrections to the linear-order 

(2) (3) 

perturbation, arising from the higher-order quantities, $0 , $q , • . • . In terms of the basic pieces of the diagrams 
shown in Fig.[ll the corrections Pab°°^{k) and Pab^°°^{k) can be diagrammatically written as the one-loop and two- 
loop diagrams, i.e., connected diagrams including one and two closed loops (e.g., see Fig. 5 in Ref. HI]), and they are 
roughly proportional to PqAq and PqAq, where Aq = k^Po{k)/{2Tr'^). The explicit expressions for the power spectra 
together with the solutions of higher-order perturbation are summarized in Appendix [Xj 

It should be noted that in the standard PT expansion, the positivity of the perturbative corrections is not guar- 
anteed. As we show later, the one- and two-loop contributions change the sign depending on the scale, and the 
absolute values of their amplitudes become comparable at lower redshift. In this respect, the standard PT has a 
poor convergence property, and the improvement of PT predictions may not be always guaranteed even including the 
higher-order corrections. 

By contrast, renormalized PT^ re-organizes the naive expansions of the standard PT by introducing the non- 
perturbative statistical quantities [s^]. In terms of these quantities, partial resummation of the naive expansion series 
is made, and the resultant convergence of the expansions is dramatically improved. In the renormalized PT, the power 
spectrum Pab{k;r]) is expressed in the form as 

Pabik;r,) = Gac{k\Ti,m)Gbd{k\Ti.m)Pcd{k;vo) + P^X''\k-ri^m) (3.4) 

with TjQ being the time at which initial condition is imposed. Here, Pcd{k; rjo) is the power spectrum given at an early 
time 770. The quantity Gab is one of the non-perturbative statistical quantities called non-linear propagator, together 
with the non-linear power spectrum. It is defined by 

(-^§^^^)^6^ik-k')Gab{\k\\v,v'); V>V', (3.5) 

where 5 stands for a functional derivative. The propagator Gab describes the infiuence of an infinitesimal disturbance 
for $a(fc'; v') on ^a{k; rj), and it coincides with the linear propagator gab in the limit fc ^ 0. Note that there is another 
non-perturbative statistical quantity called full vertex, which is the non-linear counterpart of the vertex function "fabc 



^ Strictly speaking, this statement is valid only when the universe at an early time is approximately described by the Einstein-de Sitter 
universe 

^ In this paper, we intend to make a clear distinction between the terms 'renormalized PT' and 'RPT'. While the renormalized PT 
indicates the general non-perturbative formalism developed by Ref. '325, the RPT is meant to imply the practical approximation 
method for computing the power spectrum based on the renormalized PT, which has been developed by Ref. [34l (see Appendix [Bt . 
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In the expression (j3.4p . the term P^^^ represents the corrections coming from the loop diagrams. In contrast to 

the standard PT, the loop diagrams in P^^^'' are whole irreducible, as the result of renormalization or re-organization. 
Further, each of the irreducible diagrams consists of the non-perturbative quantities of non-linear power spectrum, 
non-linear propagator and full vertex. In this respect, renormalized PT is a fully non-perturbative formulation, and 
even the expansions truncated at some levels still contain the higher-order effects of non-linear gravitational evolution. 
This is the basic reason why the convergence properties in the renormalized PT are expected to be improved. As a 
trade-off, however, a straightforward application of renormalized PT seems difhcult because of its non-perturbative 
formulation. While the term P^'~^ collects only the irreducible diagrams, it is expressed as an infinite sum of the loop 
diagrams, each of which involve the non-linear power spectrum itself. In practice, the approximation or simplification 
is needed to evaluate the expressions (j3.4p . which we will discuss in next subsection. 



B. Closure approximation 



In this subsection, taking a great advantage of the formulation of renormalized PT, we discuss how to approximately 
treat Eq. ()3.4p without losing its non-perturbative aspect as much as possible. 

In the framework of renormalized PT, the non-perturbative effects on the power spectrum are largely attributed 
to the non-linear propagator. Thus, it seems essential to give a framework to treat both the non-linear propagator 
and power spectrum on an equal footing. As it has been pointed out by Ref. [3^ . a similar kind of the renormalized 
expansion to the power spectrum (|3.4p can be made for the non-linear propagator: 

Gafc(fc|?7, v') = 9abiv, v') + r], v), (3-6) 

where the term G^^*^^ represents the mode-coupling correction, which is also made of the infinite sum of irreducible 
loop diagrams. 

In order to give a self-consistent treatment for both Eqs. (|3.4p and (|3.6p . a simple but transparent approach is to 
first (i) adopt the tree-level approximation of the full vertex function, and to (ii) apply the truncation procedure to 
the mode-coupling terms. This treatment has been frequently used in the statistical theory of turbulence in order to 
deal with the Navier-stokes equation, and is called closure approximation (43| . In the first approximation (i), the full 
vertex function is simply replaced with the linear-order one, i.e., 'jabc defined in Eq. (|2.4p . As for the truncation (ii), 
the simplest choice is to keep the one-loop renormalized diagram only, and to discard all other contributions. 

With this approximation, the mode-coupling terms in Pab and Gab are simply described by P^^^"^ — Pa^'~^'^ loop) 
and g'"^'^^ ~ ^IT^'^ loop)^ rp-j^g analytical expressions for the one-loop contributions becomes [i^ 

p^Mci-ioop)^^.^^^,^ = r dvi r dv2Gac{k\v,m)Gbd{kW,V2)^cd{k;v2,Vi), (3.7) 

Jrio Jr/o 

■n rm 



G^^^'^~^°°^\k]ri,r]') = drji dr]2 gadv^Vi) Gsb{k\T]2,v') 

f d^Q 

X 4y j^;^Jcpqiq,k ~ q) Ppriq;r]i,T]2) Gqi{\k - q\\r]i,r]2)-firs{-q,k). (3.8) 

The integrand in P^^*^'^ ioop) ^^^^^^^^^ ^j^^ function $(fc; t^i, 772), which represents the non-linear mode-coupling 
between different Fourier modes, given by 

f d^q 

$a6(fc;?7i,7?2) 2 j j^—^'^ars{q,k-q)'ybpq{q,k~q) 

X [Ppriq;m,m)Pqs{\k - q\;m,m) ^ivi - m) + Prpiq;v2,m)Psq{\k - q|; 772,771) 6(772 - ?7i)}. (3.9) 

Note that the mode-coupling function $ possesses the following symmetry: ^ab{k\r}i,ri2) — ^ba{k]rj2,rji). The 
corresponding diagrams to the integral expressions for power spectrum and non-linear propagator, i.e., Eqs. I|3.4p and 
p.6p with mode-coupling terms ((3J)) and ([31]), are shown in Fig. H 

It is worth mentioning that the integral equations ()3.4p and ()3.6p with truncated mode-coupling terms ()3.7p and 
()3.8p can be recast in the form of the integro-differential equations, and both the power spectrum and non-linear 
propagator can be computed by solving the evolution equations. This forward treatment seems especially suited 
for the full non-linear treatment of closure approximation and would be faster than directly treating the integral 



6 




FIG. 2: Diagrammatic representation of the power spectrum and non-linear propagator in closure approximation. The thick 
lines represent the full-order quantities, while the thin line indicates the linear-order one. The second terms at right-hand 
side indicate the irreducible one-loop diagrams of the mode-coupling terms, f^^*"^'^ and G^^'~^'^ loop) ^j^^ renormalized 
PT, the mode-coupling term is expressed as an infinite sum of the irreducible loop corrections. Truncating the infinite sum 
at one-loop order and adopting the tree-level approximation of the full vertex function, we obtain the closed system of power 
spectrum and propagator, as shown in the figure. 

equations. Numerica l alg orithm to solve evolution equations, together with preliminary results, is presented in details 
in Ref. ^ (see also HI). 

In the present paper, we are especially concerned with the evolution of BAOs around k < 0.4/iMpc~^, where the non- 
linearity of gravitational clustering is rather mild, and the analytical treatment even involving some approximations is 
still useful. Here, employing the Born approximation, we analytically evaluate the integral equations p.4p and (j3.7p 
[i^ . A fully numerical study on BAOs without Born approximation will be discussed in a separate paper. 

The Born approximation is the iterative approximation scheme in which the leading-order solutions are first obtained 
by replacing the quantities in the non-linear integral terms with linear-order ones. The solutions can be improved 
by repeating the iterative substitution of the leading-order solutions into the non-linear integral terms. Consider the 
time evolution of the power spectrum started from the time rjQ. For a sufficiently small value of 770, the early-time 
evolution of power spectrum is well-approximated by the linear theory. Assuming the growing-mode initial condition, 
we have 

Pabik;vo)^e^'^'>UaUi,Poik) (3.10) 

with Ua ~ (1, 1). Then, substituting Eq. (|3.10p into (|3.4p . the iterative evaluation of the the integral equations l|3.4p 
with (|3.7p by the Born approximation leads to [1^ 

Pab{k;r^) = G„(fc|r7,%)Gfc(fc|77,r;o)e2"°Po(fc) + + ' ' ' . (3-11) 

where we define Ga = Gai+Ga2- The terms P^^"^^^ and -Pq^^^"* respectively represent the leading- and next-to-leading 
order results of the Born approximation to the mode-coupling term p.7p . The explicit expressions become 

/j3 
^ laik, q; 77, ,70) h{k, q; 77, ,70) e'^«Po{q) Po{\k - q|), (3.12) 

X e^"" Po(|fe - p\)Po{q)Poi\p - q\). (3.13) 

The kernels la and Jq are respectively given by 

Ia{k, q;r], 110) ^ f dr]' Gai{k\T],i]') -yirsiq, k - q) Gr{qW ,Vo) Gsi\k - q\\T]' ,t]o), (3.14) 
J no 

Ja{k, p, q;r], r]Q) = i dr]i / d'q2Gai{k\r],r]i) jirsip, k - p) Grcip\vi,V2) 

Jijo Jvo 

X Icpqiq^P - q)Gp{q\r]2,'no)Gqi\p ~ q\\m,Vo)Gsi\k - pll-qi^Tjo)- (3.15) 




FIG. 3: Diagrammatic representation for the perturbative treatment of the power spectrum with the Born approximation, i.e., 
Eq. dMll- 



The diagram corresponding to the above expressions is shown in Fig. [31 Note that in deriving the expression 
(|3.1ip . we do not expand the propagators Gab and their non-perturbative properties stiU hold. In order to evaluate 
Eq. (|3.1ip . we use the analytic solution of Gab derived in Ref. [43*1, where the non- linear propagator was constructed 
approximately by matching the asymptotic behaviors at low- and high-k modes, based on Eqs. (j3.6p with p.8p . The 
resultant analytic solution behaves like Gab — > Sab Ji{2x)/x at fc — > oo, where the quantity Ji is the Bessel function 
with its argument x — k(7v{e^ ~ e'' ), and the velocity dispersion is approximately described by the linear theory, 
i.e., ~ lin — J dq Piin{q; z) / {dir'^) . Note that the final results of the power spectrum are a little bit sensitive 
to the high-fc behavior of the propagator, and a naive application of the approximate solution leads to a slight shift 
in the amplitude of power spectrum. While this is not serious at all for the leading-order calculation, it amounts 

(MC2) 

a percent-level shift when we consider the higher-order correction, P^^ . As discussed by Crocce & Scoccimarro 
(2008), one possible reason for this may be a small contribution from the sub- leading corrections in the propagator. 
In order to remedy the effect of small corrections, we follow the method proposed by Ref. 0|. We define 



a{z) = 



Jg™ dqPni{q;z) 



-, 1/2 



(3.16) 



where P^i means the non-linear matter power spectrum. Then, the sub-leading correction can be corrected by simply 
multiplying the factor a by a^, i.e., — > a{z)(7v Note that this treatment is only applied to the propagator in 
the lowest-order term in Eq. (|3.11|) . which most sensitively affects the power spectrum amplitude on small scales. 
For simplicity, we use halof it 50] to compute P^i and adopt the cutoff wavenumber, fcmax = k^, where fco- is the 
non-linear scale defined by Ref. i50l] . 

In the rest of this paper, we present the results for the analytic treatment based on the expression (|3.1ip . In 

computing the mode-coupling terms Pa^'~^^'^ and P^^'~^^\ we must first evaluate the functions la and Jq for a given 
set of arguments, which involve the one- and two-dimensional integrals over time rj. We use the Gaussian quadrature 
for these time integrations. As for the momentum integrals in the mode-coupling terms, thanks to the symmetry 
of the functions la and Ja, the multi-dimensional integrals in P^^'^^'' and P^^*^^' can be reduced to the two- and 
four-dimensional integrals, respectively. We use the Gaussian quadratures for the momcntmn integral in P^^'^^''. The 

four-dimensional momentum integration in the mode-coupling term is performed with Monte Carlo technique 

of quasi-random sampling using the library, Cubans Ij"^. 

Finally, we note that the formulation and analytic treatment presented here have several distinctions and similarities 
to the other non-perturbative calculations proposed recently. In Appendix [Bl we compare the present work with a 
subset of these treatments, and discuss how the approach developed here is complementary to or expands on these 
studies. 



IV. IMPROVED PT VS. NUMERICAL SIMULATIONS 



In this section, particularly focusing on the BAOs, we compare the improved PT predictions from the analytic 
treatment of closure approximation with results of N-body simulations. 



^ |http : //mw . f eynarts . de/cubaZ 
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A. N-body simulations 

We use a publicly available cosmological N-body code, Gagdet2 [s^]- We ran two sets of simulations, wmap3 and 
wmap5, in which we adopt the standard Lambda CDM model with cosmological parameters determined from the 
WMAP3 and WMAP5, respectively The wmap3 run is basically the same N-body run as described in Ref. [45| , 

and a quantitative comparison between the leading-order results of improved PT and simulations has been previously 
made. We basically use the results of wmap3 run to check the consistency of the present calculations with the previous 
work. The wmapS run is also helpful to cross-check the convergence properties in the new simulation, wmap5, which 
increase the number of realizations to 30. Table [T] summarizes the parameters used in the simulations. The initial 
conditions were created with the 2LPT code [53] at initial redshift Zi^i = 31, based on the linear transfer function 
calculated from CAMB [54]. The number of meshes used in the particle-mesh computation is 1,024'^. We adopt a 
softening length of O.lft-^^Mpc for tree forces. 

We store three output redshifts for wmapS run, whereas we select four output redshifts for wmap5 run; z = 3, 1, 
and (wmap3) : z = 3, 2, 1, and 0.5 (wmap5). Using these outputs, we compute the power spectrum and two-point 
correlation function in both real and redshift spaces. 

The calculation of the matter power spectrum adopted here is basically the same treatment as in Ref. |45j . The 
standard method to compute the power spectrum is to square the Fourier transform of the density field and to take 
an average over realizations and Fourier modes. This is given by 

^(fc") = WN^n E E \s--'\kf ; K^j^ y: ' (4-1) 

" ™=1 fc»i"< I fc| " *:™<|fc|<r»='- 

where and TV^" are the number of Fourier modes in the n-th wavenumber bin and the number of realizations, 
and /c™'" and k^'^^ are the minimum and the maximum wavenumber of the n-th bin, respectively. The quantity 
(5™~"^(fc) means the density field in Fourier space obtained from the m-th realization data. We use the Cloud-in-Cells 
interpolation for the density assignment of particles onto a 1,024^ mesh, and correct the window function. Note 
that the power spectra measured from the standard treatment above suffer from the effect of finite-mode sampling 
discussed by Ref. ^5^. The resultant power spectrum deviates from the prediction for the ideal ensemble average, and 
exhibits the anomalous growth of power spectrum amplitude on large scales. In order to reduce the effect of finite 
mode sampling at fc < 0.1/iMpc~^, we multiply the measured power spectrum by the ratio, P^"^ {k)/ P\in{k), where 
the quantity P^^{k) is calculated from the perturbation theory up to the third-order in density field, and Piin{k) is 
the input linear power spectrum extrapolated to a given output redshift. Note that in computing P^^{k), we use 
the Gaussian-sampled density field used to generate the initial condition of each N-body run. With this treatment, 
the individual random nature of each N-body run is weakened, and the errors associated with anomalous growth is 
reduced"* . 

For the estimation of two-point correlation function, we adopt the grid-based calculation using the Fast Fourier 
Transformation (FFT). In this treatment, similar to the power spectrum analysis, we first compute the square of the 
density field on each grid of Fourier space. Then, applying the inverse Fourier transformation, we take the average 
over realization and distance, and obtain the two-point correlation function. Schematically, this is expressed as 

^(-") = ]^;]V^ E E [l<J"-^fe)P;r], (4.2) 

" m=l rn"n<|r|<rn^'»'' 

where the operation FFT stands for the inverse FFT of the squared density field on each grid. Note here that r„ 
is simply chosen at the center of the n-th radial bin, i.e., r„ = {ryain + ''max)/2. 

Eq. (14. 2p usually suffers from the ambiguity of the zero-point normalization in the amplitude of two-point correlation 
function, because of the lack of the low-fc powers due to the finite boxsize of the simulations. With the 1, 024'^ grids 
and the boxsize of Lbox — l/i~^Gpc, however, we can safely evaluate the two-point correlation function around the 



* In Ref. [45| . the correction to the effect of finite-mode sampling has been applied to the real-space power spectra. Itere, we extend it 
to compute the redshift-space power spectrum by simply replacing the ratio P^"^ {k) / Puj^lk) with that in redshift space. To be precise, 
we compute the multipole moments of the redshift-space power spectrum, and the ratio, (A:)/P^ jj^(fc), is multiplied for each 

multipole spectrum (see Sec. IIV C 11 1. 
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TABLE I: Parameters of N-body simulations 



Name 




# of particles 




^ of runs 








h 


ris 




wmap3 


1000/i-^Mpc 


512^ 


31 


4 


0.234 


0.766 


0.175 


0.734 


0.961 


0.76 


wmap5 


1000/i"^Mpc 


512^ 


31 


30 


0.279 


0.721 


0.165 


0.701 


0.96 


0.817 




10-2 10-1 iQ-z 10-1 

k [h Mpc-'] k [h Mpc-'] 



FIG. 4: Convergence properties of standard PT (left) and improved PT (right) expansions in the matter power spectrum. 
In each panel, the higher-order contributions to the total power spectrum labeled as Pni is separately plotted. In left panel, 
one- loop and two-loop corrections in the standard PT, p^_^^°°P and P^\'°°^, are plotted, while in right panel, the mode-coupling 
corrections p^^'-"^^ and p^'^*-'^^ in the improved PT given at Eqs. (|3.12|) and p.l3|) respectively shown (labeled as MCI and 
MC2), together with the first term in Eq. (f3TTT|l (labeled as G^Po ). Note that dashed lines indicate the negative values. 

baryon acoustic peak. Comparison between different computational methods, together with convergence check of this 
method, is presented in Appendix [Cl 

Finally, similar to the estimation of power spectrum, the finite-mode sampling also affects the calculation of the 
two-point correlation function. We thus correct it by subtracting and adding the extrapolated linear density field as, 
S.(r) — ^iin('') + Ciin(f), where is the correlation function estimated from the Gaussian density field, and ^lin is the 
linear theory prediction of two-point correlation function. 



B. Results in real space 



1. Power spectrum 



Before addressing a quantitative comparison between N-body simulation and improved PT, we first discuss the 
convergence properties of the improved PT, and consider how well the calculation based on the improved PT does 
improve the prediction compared to the standard PT. 

Fig.Hlplots the overall behaviors of the non-linear power spectrum of density fluctuation, P{k; z) = Pii(fc; z), given 
at z = 0, adopting the wmap3 cosmological parameters. In left panel, the results of standard PT are shown, and the 
contributions to the total power spectrum up to the two-loop diagrams are separately plotted. On the other hand, 
right panel shows the results of improved PT. We plot the contributions up to the second-order Born approximation 
labeled as MCI and MC2. 

In Fig. 01 there are clear distinctions between standard and improved PTs. While the loop corrections in standard 
PT change their signs depending on the scales and exhibit an oscillatory feature, the corrections coming from the Born 
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FIG. 5: Ratios of power spectrum to smoothed reference spectrum, P(fc)/Pno-wiggic(fc), given at redshifts z — 3(top), l(middle) 
and O(bottom). Cosmological parameters used in the wmap3 simulations are adopted to compute the power spectrum from 
standard PT and improved PT, and the resuhs are compared with N-body simulations (symbols with error-bars). The reference 
spectrum Pno-wiggic (fc) is calculated from the no- wiggle formula of the linear transfer function in Ref. [2^. In each panel, dotted, 
dashed and solid lines represent the linear, standard PT and improved PT results, respectively. In left panel, leading-order 
results of standard PT and improved PT are shown, while in right panel, the results including the higher-order corrections are 
plotted. 



approximation in the improved PT are all positive and mostly the smooth function of k. Further, the higher-order 
corrections in the improved PT have a remarkable scale-dependent property compared to those in the standard PT; 
their contributions are well-localized around some characteristic wavenumbers, and they are shifted to the higher k 
modes as increasing the order of PT. These trends clearly indicate that the improved PT with closure approximation 
has a better convergence property. Qualitative behaviors of the higher-order corrections quite resemble the predictions 
of RPT by Crocce & Scoccimarro (2008) [sl. 

Now, let us focus on the behavior of B AOs, and discuss how the convergence properties seen in Fig. 2] affect the pre- 
dictions of BAO features. In Fig.[5l adopting the wmap3 cosmological parameters, we plot the ratio, P(fc)/Pno-wiggio(fc), 
where the function Pno-wiggic(fc) is the linear power spectrum from the smooth transfer function neglecting the BAO 
feature in Ref. [23| . In left panel, N-body simulations are compared with the leading-order results of PT predictions, 
i.e., standard PT including the one-loop correction (dashed), and improved PT with first-order Born correction (solid). 
Apart from the wiggle structure, the amplitude of standard PT prediction monotonically increases with wavenumber 
k, and tends to overestimate the results of N-body simulations. On the other hand, the amplitude of improved PT 
prediction rapidly falls off at a certain wavenumber, and the deviation from N-body results becomes significant. How- 
ever, a closer look at the behavior on large scales reveals that improved PT prediction gives a better agreement with 
simulation. The results are indeed consistent with the previous findings in Ref. [45j . The situation becomes more im- 
pressive when we add the next-to- leading order corrections. As shown in right panel, the improved PT gets the power 
on smaller scales, and reproduces the N-body results in a wider range of wavenumber. By contrast, the prediction of 
standard PT depicted as dashed lines seems a little bit subtle. Compared to the one-loop results, the amplitudes of 
the standard PT prediction including the two-loop correction are slightly reduced, and the agreement with N-body 
simulation seems apparently improved a bit at higher redshift. At lower redshift z = 0, however, the correction coming 
from the two- loop order becomes significant, and the prediction eventually underestimates the simulation. The reason 
for these behaviors basically comes from the competition between positive and negative contributions of the one-loop 
and two-loop corrections, respectively (see left panel of Fig . [3]) . These are consistent with those findings in Ref. [40| 
(see Fig. 1 of their paper). 

In Fig. [51 to clarify the range of agreement in more quantitative ways, we plot the fractional difference divided by 
the smoothed reference spectra, [PN-body(fc) ~ ^'pT(fc)]/^'no-wiggio, where the quantity PpT(fc) implies the standard 
and improved PT predictions in left and right panels, respectively. Here, the vertical arrows represent the maximum 
wavenumber ki%, below which the leading-order predictions of standard or improved PT reproduce the N-body results 
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FIG. 6: Difference between N-body and PT results divided by the reference spectrum, [PN-body(fc) — PpT{k)]/ Pno-wiggia{k) . Left 
panel shows the results for standard PT up to the two-loop order. Right panel presents the case of improved PT including the 
corrections up to the second-order Born approximation of the mode-coupling term. In both panels, vertical arrows represent the 
wavenumbers ki% of standard and improved PT (from left to right), below which the leading-order PT predictions reproduce 
the N-body simulations well within 1% accuracy (see text in details). 



quite well within the 1% accuracy. According to Nishimichi et al. [4^, this has been determined by the detailed 
comparison between models and simulations, and is empirically characterized by solving the following equation: 



A-2 



(4.3) 



with C ~ 0.18 for the one-loop standard PT, and C — 0.35 for the improved PT up to the first-order Born correction. 

Comparing these convergence regimes of the leading-order calculation with results of fractional differences. Fig. [H] 
shows that the inclusion of higher-order terms does not always improve the prediction in the standard PT treatment. 
By contrast, the improved PT calculation does improve the predictions, and the range of agreement between N-body 
simulations and the predictions becomes wider. 

In Fig. [3 we plot the results for the wmap5 simulations, which have relatively large value of erg compared to the 
wmap3 run (see Table H]). Left and right panels respectively plot the ratio of power spectrum amplitude and the 
fractional difference between N-body results and improved PT predictions. With the 30 runs of N-body simulations, 
the errors in the power spectrum amplitude are greatly reduced, and it is clearly shown that the predictions of 
improved PT including the higher-order corrections almost coincide with the N-body results beyond the convergence 
regime of the leading-order calculations (indicated by vertical arrows), and achieve a sub-percent accuracy. From this 
plot, the maximum wavenumber ki% at each redshift can be estimated by comparing the predictions with N-body 
results as ki% = 0.20/iMpc-^(z=0.5), 0.23/iMpc-i(z=l), 0.33/iMpc-i(z=2) and 0.47/iMpc-i(z=3). These values 
roughly match those determined from the criterion (|4.3p with the constant C = 0.70. 

Although we did not store the z = data of wmap5 run to compare with analytic prediction, Eq. (|4.3p using this 
constant value implies that the maximum wavenumber for improved PT becomes ki% = 0.15/iMpc~^, which contrasts 
with the one for the one-loop prediction of standard PT, ki% = 0.09/iMpc^^. Thus, the improved PT including up to 
the second-order Born approximation is expected to be still accurate at z = 0, and it can cover the major part of the 
BAOs. A more detailed comparison at low redshift including other analytic prescriptions can be found in Ref. [46j . 
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FIG. 7: Comparison between N-body results and improved PT predictions in the case adopting wmapS cosmological parameters. 
From top to bottom, the results at z = 3, 2, 1 and 0.5 are shown. The improved PT predictions plotted here include the 
corrections up to the second-order Born approximation of the mode-coupling term, p^'^^ Lgf^. j-atio of power spectrum to 
the smoothed reference spectra, P(fc)/Pno-wiggic(fc). Solid and dotted lines are improved PT and linear theory predictions, 
respectively. Right: difference between N-body and improved PT results normalized by the no-wiggle formula, [PN-body(fc) — 
PpT(fc)]/Pno-wiggio(fc). In each panel, vertical arrows represent the wavenumber ki% for the leading-order predictions of standard 
and improved PT (from left to right). 



2. Correlation function 



Having confirmed the excellent properties of the improved PT, we turn to focus on the baryon acoustic peak in the 
two-point correlation function. The two-point correlation function can be computed from the power spectrum as 

^, , f dkP ^ ,,,sin(fcr) ,^ ^, 

^ir)-J ^Pnik)^. (4.4) 

Top panel of Fig. [5] shows the two-point correlation functions around the baryon acoustic peak at different redshifts 
z = 0.5, 1,2 and 3 (from top to bottom) in the case adopting the wmap5 cosmological parameters. Also, lower panel 
plots the fractional differences between N-body and improved PT results, i.e., [•fN-body (?") — Cpt(?')]/Cpt('')- 

After the correction of finite-mode sampling, the error-bars in N-body simulations are greatly reduced, and the 
deviation of the N-body results from linear theory predictions (depicted as dotted lines) is clearly seen. As decreasing 
the redshift, the baryon acoustic peaks become smeared and the position of the peak are slightly shifted to a smaller 
scale. These trends can be accurately described by the leading-order calculation of improved PT, and the agreement 
between N-body results and the predictions is excellent. The fractional error in amplitude is well within a few percent, 
except for a large separation beyond the location of baryon acoustic peak, where the accuracy of N-body results tends 
to be worsen due to the limited simulation boxsize. Note that the corrections coming from the higher-order Born 
approximation do not alter the behaviors at r > 30/i~^Mpc, and their amplitudes are negligibly small compared to the 
error-bars of N-body simulations. Thus the leading-order prediction seems robust for describing the baryon acoustic 
peak. 

It has been recently suggested by several authors that the smearing effect on baryon acoustic peak is mostly 
attributed to the random motion of mass distribution [Hll, and it is approximately described by the convolution of 
the Gaussian smoothing function (e.g., [s^, \E^). In the language of improved PT, this effect corresponds to the 
disappearance of the memory of initial condition, which is encoded in the non-linear propagator. Strictly speaking, 
the asymptotic behavior of the non-linear propagator is not a Gaussian form in closure approximation, although the 
damping behavior manifestly exhibits in the approximate solution of non-linear propagator. Hence, the prediction for 
the two-point correlation function seems robust against the high-fc behavior of the non-linear propagator. 

Finally, it should be noted that the standard PT prediction fails to converge the integral in Eq. (|4.4p . because of 
the high-fc behavior of the power spectrum. This is true even when including the higher-order correction of two-loop 
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FIG. 8; Top: Two-point correlation functions in real space adopting the wmap5 cosmological parameters. The solid lines 
represent the leading-order predictions of improved PT, while the dotted lines show the linear theory results. Bottom: Fractional 
differences between N-body and improved PT results, [^N-body(7") — Cpt(?")]/^pt(?'). In both panels, the symbols with error-bars 
indicate the N-body results averaged over the 30 realizations in which the effect of finite-mode sampling is corrected: z = 0.5 
(open stars), 1 (open squares), 2 (filled triangles), and 3 (crosses). 

order. Thus, the successful results of improved PT prediction may be regarded as an outcome of non-perturbative 
property. 

C. Results in redshift space 

In practical observation with galaxy redshift surveys, the observed galaxy distribution is inevitably distorted due 
to the peculiar velocity of each galaxy. The so-called redshift-space distortion is known to alter the shape of the 
power spectrum in two different ways (e.g., [58|). One is the apparent enhancement of the clustering signal called 
Kaiser effect [H^, which originates from the bulk motion of mass distribution falling into the massive halos. Another 
important effect is the finger-of-God (FoG) effect, which effectively suppresses the power spectrum amplitude on small 
scales by the virialized random motion of the mass residing at a halos. 

Although a rigorous non-perturbative treatment of the redshift-space distortion is difficult, these two effects has 
been phenomenologically modeled as (e.g., [60l. [6ll. [g^. [b^ ) 

Pfs^fc, ^i) = {l + fi^ ff Pll(fc) i^FoG(fc m), (4.5) 

where (x is the cosine of the angle between the line-of-sight direction and the Fourier mode fc, and / is the logarithmic 
derivative of linear growth factor, defined as / = d\nD / dlna. The function DfoG represents the damping function 
which mimics the FoG effect, and it asymptotically approaches unity in the fc ^ limit, where the linear-theory 
formula by Kaiser is recovered. 

Recently, Scoccimarro [g^l proposed an improved version of the model (14.51) to properly take account of the non- 
linear evolution of density and velocity fields on the Kaiser effect (see also [65l [66|). This is expressed as 




p(S)(fc,^) = [Pn(fc) + 2//i2pi2(fc) + /V-P22(fc)] exp{-(/MA:av)'}. 



(4.6) 
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FIG. 9; Ratio of power spectra to smoothed reference spectra in redshift space, (^)/^ino-wiggio(^)' fro™ the wmap5 simu- 

(s) 

lations. The reference spectrum no-wiggie calculated from the no-wiggle approximation of the linear transfer function, and 
the linear theory of the Kaiser effect is taken into account. Left panel shows the monopole power spectra = 0), and the right 
panel shows the quadrupole spectra {i = 2). Solid and dashed lines represent the results from the improved PT adopting the 
model of redshift-space distortion (|4.6|) . To plot the results, the linear theory was used to compute ctv in dashed lines, while in 
solid lines, ctv was determined by fitting the predictions to the N-body simulations. In each panel, vertical arrow indicates the 
maximum wavenumber ki% for improved PT prediction including up to the second-order Born approximation, which has been 
estimated from Fig. [7] (see Sec. IIV BTI for definition of ki%). 



Here, the quantity CTv is the one-dimensional velocity dispersion given by 

In what follows, we adopt the model (|4.6p to calculate the redshift-space power spectrum. Although this model is 
still phenomenological and may not be regarded as the best one, a comparison between the model predictions and 
N-body simulations shows that the prediction based on the model (|4.6p gives a better result. Taking Eq. (|4.6p as 
a canonical model of the redshift-space distortion, we will investigate the extent to which the model (|4.6[) faithfully 
reproduces the N-body results well, and discuss how to improve the model prescription. 



1. Power spectrum 



For a quantitative comparison of model prediction with N-body simulation, we compute the multipole moments of 
the two-dimensional power spectrum P'^^^(fc, /i): 

Pf )(fc) ^ ^i±l f d^,p^^\k,n)v,{^,), (4.8) 



with Vi being the Legendre polynomials. 

Substituting the model (|4.6p into the above, the monopole, quadrupole and hexadecapole contribution to the 
redshift-space power spectrum are analytically expressed as 

PPik) = po(fc), (4.9) 
PPik) = ^{3pi(fc)-Po(fc)}, (4.10) 

Pfik) = l{35p2{k)-30pi(k)+3po{k)}. (4.11) 



15 



500 
400 
'to 300 

e 

' 200 
100 


12 3 4 

z 

FIG. 10: Redshift evolution of velocity dispersion CTv. While the solid lines represent the linear theory prediction, the open 
squares indicate the results obtained by fitting the model (|4.6|) to the monopole and quadrupole spectra of N-body simulations 
(see Fig. O. 

where the function p„(fc) is defined by 
Pn{k) = i 

The quantity 7(n, k) is the incomplete gamma function of the first kind 

7(n,K) ^ f dt r-^e"* (4.13) 
Jo 

with its argument k = {k f (Tv)^. 

Fig. [HI shows the monopole (left) and quadrupole (right) moments of the redshift-space power spectra at different 
redshifts, obtained from the wmap5 simulations. We do not plot here the hexadecapole contributions, because the 
power spectrum estimated from the N-body simulations is still noisy even with the 30 realizations. In each panel of 
Fig. [H the dashed lines indicate the improved FT predictions based on the model (j4.6p , where the corrections up to 
the second-order Born approximations are included. Note that the velocity dispersion is computed from the linear 
theory. Clearly, the predictions all underestimate the N-body results, and the agreement between predictions and 
N-body simulations is restricted to a quite narrow range on large scales. As a reference, we also show the maximum 
wavenumber fc]^% of the improved FT prediction (vertical arrows), in which we include the corrections up to the 
second-order Born approximation in real space (see Fig. [7]). 

The reason why the prediction generically underestimates the N-body simulations would be partly attributed to 
the calculation of the velocity dispersion using the linear theory. It has been advocated by several authors that 
the suppression of power spectrum by FoG effect is originated from the non-linear structure of virialized halos, and 
thereby the linear theory estimation of CTv may be inappropriate. In this respect, we admittedly regard CTv as an 
uncontrollable parameter, which should be determined by fitting the predictions to N-body results. 

The solid lines in each panel of Fig. [9] show the results of redshift-space spectra adopting the fitted values of 
CTv In estimating a^, both the monopole and quadrupole spectra were fitted to the N-body results in the range of 
< fc < fci%. Fig. 1101 summarizes the fitted results of Uv, which significantly deviate from the linear theory prediction 
at higher redshifts. 

Then, apparently, overall agreement between prediction and simulation becomes fairly improved, although as a 
trade-off, small discrepancy manifests at low-k mode, where the N-body results rather agree well with the prediction 
adopting Uv calculated from linear theory. In Fig. [Til left and right panels respectively plot the fractional differences 
of the monopole and quadrupole moments between the model predictions and N-body simulations. Except for the 
narrow range of low-k modes, a percent-level agreement is almost achieved for the monopole power spectrum. This 
is true at least within the convergence regime calibrated in real space (see vertical arrows in Fig. Ilip . However, the 




7(n+l/2,«) 



Pll(fc) 



7(n + 3/2,K) 



fPi2{k) 



7(n-|-5/2,K) 

Zn+5/2 



(4.12) 
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FIG. 11: Difference between N-body and PT results divided by the reference spectrum in redshift space, i.e., [^Y^N-bodyC^) ~ 
^ii=T(^)l/^f no-wiggicC^)- The left and right panels respectively represent the results from monopole and quadrupole power 
spectra. Note that the improved PT predictions are computed based on the model (|4.6|l adopting the fitted value of CTv. 
For comparison, the statistical errors limited by the cosmic variance of the survey volumes roughly corresponding to those of 
WFMOS-like survevj28| and B0SS[2^ are shown as shaded regions in panels of 2 = 3, 2 = 1 and z = 0.5, assuming respectively 
the survey volumes of F = lh~^Gpc^ , 4/i~^Gpc^ and 4.5/i~^Gpc^. Note that in each panel, vertical arrow indicates the 
maximum wavenumber fci% determined from Fig. [7| by comparison between N-body and improved PT results. 



fractional error of the quadrupole power spectrum still exhibits a little bit large discrepancy, signaling the fact that 
the model (|4.6p misses something important for higher-multipole moment of redshift-space distortion. 

To see the significance of this deviation in practice, in Fig. llli the expected l-cr errors limited by the cosmic variance, 

(s) 

AP^ (k), are shown, depicted as the shaded region. Here, we specifically consider the ground-based BAO surveys 
like WFMOS surveyfH and BOSSjlil, assuming the survey volumes of = 1 h-^Gpc^ at z = 3 and Ah-^Gpc^ at 
z = 1 for WFMOS survey, and V = 4.5 h~^Gpc^ at z = 0.5 for BOSS^. Based on the approximation that the density 
field is well-described by a Gaussian random field, the cosmic- variance error can be estimated as 

[APPikr^^alAk), (4.14) 

where the quantity N^. is the number of Fourier modes within a given bin at fc, and is given by = 
47rfc2A/fc/(27r/ibox)V2 ^VP Ak/{2TTf. The function ap^i is 

= i^^±ll! l^/i [p'^^\k,i,)VM]\ (4.15) 

The expression (|4.14l) with (|4.15p is a generalization of the cosmic- variance error in real space (e.g., [47l. WA. Issl, [69j ) 
to the multipole moments in redshift space. Note that the error AP^ (fc) depends on the bin width Afc, for which 
we simply adopt the same bin size as used in the power spectrum analysis of N-body data. The analytic estimate 
of ApP based on Eq. (|iTil) is roughly consistent with the statistical errors estimated from the N-body data of 30 
realizations. 

Comparison between the cosmic-variance errors and fractional differences shows that the discrepancy seen in the 
quadrupole power spectrum is definitely large, and it eventually exceeds the statistical error at large k modes. Since 



^ Strictly speaking, BOSS project is a part of Sloan Digital Sky Survey III, aiming at precisely measuring the cosmological distance and 
expansion rate at 2; = 0.35, 0.6 and z = 2.5. Here, we only consider the low-2 measurement with survey depth 0.2 ^ 2 ^ 0.8. 
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this has happened inside the vahd range of the improved PT calibrated in real space (indicated as vertical arrows) , we 
conclude that the current model prediction with (|4.6p is insufficient to describe the higher- multipole moments of B AOs, 
and a more elaborate work on the models of redshift-space distortion is needed for upcoming BAO measurement. 



2,. Correlation Junction 

Finally, we discuss the correlation functions in rcdshift space. Similar to the power spectrum, we apply the multipole 
expansion to the anisotropic two-point correlation function as 

- J P(^nA,M)e''^" = E (^)^H-) (4.16) 

£:cvcn 

with v = The multipole moment of the correlation function, Q , is directly related to the Fourier counterpart, 

P^^^ through 

^f\^) = ^' j^-^Pf\k)u{ks). (4.17) 

Fig. [T2l shows the monopole (left), quadrupole (middle) and hexadecapole (right) moments of correlation function. 
In each panel, the N-body results are compared with the predictions from linear theory (dotted) and the leading- 
order calculation of improved PT (solid) adopting the model (|4.6p with linear theory prediction of ctv Note that the 
predictions of improved PT are hardly changed by including the higher-order corrections and/or using the fitted value 
of (7v, at least around the baryon acoustic peak, and the systematic differences between including and ignoring the 
corrections are well within the error-bars of N-body simulations. 

As anticipated from the results in real space, the baryon acoustic peaks in the monopole moment tend to be 
smeared as decreasing redshift, but the effect seems little bit stronger than those in real space. This is due to the 
additional effect coming from the redshift distortion. Although no prominent signal of the BAOs exists in the higher- 
multipole moments, the same tendencies can be seen in the quadrupole and hexadecapole moments. The improved 
PT calculations are broadly consistent with N-body results, but small discrepancies manifest around the baryon 
acoustic peak and trough. Lower panels of Fig. [T^] showing the fractional differences imply that these are at most 
5% effect in amplitude, except for the hexadecapole case with large error-bars of simulation. It is interesting to note 
that no noticeable redshift dependence appears in the fractional differences, indicating that the discrepancies may be 
attributed to the model of redshift-space distortion. Furthermore, it turns out that these are well within the cosmic- 
variance errors of the ground-based BAO measurement, indicated as shaded region. Assuming that the underlying 
density field is well described by a Gaussian random field, the cosmic variance for the multipole correlation functions 
^f'^ can be written as (see [13, [zO, [7l| for cosmic-variance errors in real space) 

Aef^(^)]' = I / '^{3dks)apAk)? (4.18) 



with up^i given by Eq. (|4.15p . Note that the analytic estimation of cosmic- variance errors A^^^-* shown in Fig. [T2l 
reproduce the N-body results quite well. 

Hence, compared to the power spectrum in redshift space, the correlation functions obtained from the N-body 
simulation and analytic calculation can have a better agreement. Presumably, this is because the acoustic peak 
structure in the correlation function is mostly attributed to the low-k behavior of the BAOs, and the power spectrum 
at low-k modes is accurately described by the model (|4.6p with the linear theory prediction of CTv In other words, 
the baryon acoustic peak would be robust against the non- linear effects at high-k modes (see also [s^ H^, [T^I ) • This 
implies that even the prediction at the current level is sufficient to characterize the acoustic peak in the correlation 
function, and it can be used as an accurate theoretical template for future precision BAO measurement. 

Note, however, that the measured amplitudes of the two-point correlation function arc strongly correlated between 
different scales. In practice, not only the diagonal component but also the off-diagonal components of the covariance 
of the correlation function must be considered for a reliable estimation of cosmological distance, and a more careful 
study is needed. 



V. DISCUSSION AND CONCLUSIONS 



In this paper, we have presented the improved PT calculations of the matter power spectrum and two-point 
correlation function in real and redshift spaces. Based on the closure approximation of the renormalized PT treatment. 
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FIG. 12: Top: Correlation function in redshift space. Left, middle and right panel respectively shows the monopole, quadrupole 
and hexadecapole contributions to the anisotropic correlation function The solid and dotted lines are the predictions from 
the improved PT based on the model (|4.6p and linear theory, respectively. Note that only the leading-order Born approximation 
to the the mode coupling term is included in the improved PT. z = 0.5 (red); z — 1 (magenta); z = 2 (cyan); z — 3 (green). 
For comparison, the statistical errors limited by the cosmic variance of the survey volumes V = l/i~^Gpc^, 4/i~^Gpc^ and 
4.5/i-^Gpc^ are estimated from Eq. (|4.18|l . and are depicted as shaded regions around the N-body results at z — 3, z = 1 
and z = 0.5, respectively. The cosmic-variance error for hexadecapole is not shown here because of the large scatter. Bottom: 
Fractional differences of the results between N-body simulations and improved PT predictions, [^N-body (s) — Cpt(s)]/Cpt(s) for 
different redshifts ai z = 0.5 (open stars), z = 1 (open squares), 2 = 2 (filled triangles), and z — 3 (crosses). 



a closed set of the non-perturbative expressions for power spectrum and propagator is obtained. The resultant 
expression includes the effect of resummation for a class of loop diagrams at infinite order, and thereby the convergence 
of higher-order contributions is expected to be improved. Employing the Born approximation, we have analytically 
calculated the non-linear power spectrum, and compared the convergence properties of improved PT with those of 
standard PT by explicitly computing the higher-order corrections. 

We have also made a detailed comparison between the improved PT result and N-body simulations. With a large 
boxsize and many realization data of N-body simulations, the statistical errors of two-point statistics are greatly 
reduced by the correction of the effect of finite-mode sampling, and this enables us to investigate the convergence 
check of numerical and analytic calculations at a percent level. Then, specifically focusing on the behaviors of BAOs, 
the power spectrum and two-point correlation functions are calculated in both real and redshift spaces. In redshift 
space, the effect of redshift-space distortion which changes the clustering pattern of mass distribution should be 
incorporated into the improved PT predictions. In this paper, adopting the model proposed by Ref. ^64(1 (Eq. (|4.6p ), 
we have quantified the extent to which the current model description faithfully reproduces the N-body results, and 
clarified the key ingredients toward an improved prescription of redshift-space distortion. 

Our important findings are summarized as follows: 

• The improved PT expansion based on the Born approximation has better convergence properties, in marked 
contrast with the standard PT expansion. The corrections coming from the mode-coupling term are well- 
localized positive functions of wavenumbers, and their contributions tend to be shifted to a higher k region 
as increasing the order of perturbation. Thus, the inclusion of higher-order corrections stably improves the 
prediction, and the range of agreement with N-body results becomes wider in wavenumber. 

• In real-space power spectrum, the improved PT prediction including up to the second-order Born correction 
seems essential for modeling BAO precisely. We estimated the maximum wavenumber below which the 
results of both the N-body simulation and improved PT calculation converge well within the 1% accuracy. The 
resultant value of fci% can be summarized as Eq. (j4.3p with the constant value C = 0.7, which provides a way to 
estimate in a cosmology independent manner. On the other hand, if we consider the two-point correlation 
function in real space, the leading-order calculation turns out to be sufficiently accurate, and no higher-order 
correction is needed to describe the non- linear evolution of baryon acoustic peak seen in the N-body simulations. 

• Modeling redshift-space power spectrum with Eq. (14. 6p gives a broadly consistent result with N-body simulations, 
if we regard the velocity dispersion as a fitting parameter. However, discrepancy between improved PT 
predictions and N-body results has appeared in the quadrupole power spectrum, and it becomes larger than the 
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statistical errors limited by the cosmic variance of the survey volume V ~ a few h~^Gpc^. This is true even in 
the valid range of improved PT, k < On the other hand, while a small descrepancy has been also found in 
the two-point correlation, it turns out that the discrepancy is well within the cosmic-variance error, and even the 
leading-order prediction using the linear theory estimate of CTv can be used as an accurate theoretical template 
for future ground-based BAO measurement. 

The recently proposed techniques to deal with the non-linear gravitational clustering, including the present treat- 
ment, have been greatly developed, and they would be a promising cosmological tool to precisely model the shape and 
amplitude of the power spectrum and/or the correlation functions in an accuracy of sub-percent level. Combining the 
model of redshift-space distortion, we are now able to discuss the non-linear clustering in redshift space. Although 
the present paper is especially concerned with the analytical work, we note that the non-perturbative formulation 
with closure approximation is suited for forward treatment in time [48| . in which all orders of Born approximation 
can be fully incorporated into the predictions by numerically solving the evolution equations. This approach would 
be particularly useful to study the non-linear matter power spectrum in general cosmological models, including the 
modified theory of gravity ^4S^ . 

Finally, in practical application to the precision BAO measurements, there are several remaining issues to be 
addressed in the future work. The improvement of the model of redshift-space distortion is, of course, a very important 
and urgent task. The effect of galaxy biasing is also one of the key ingredients for modeling accurate theoretical 
template, and several attempts to take account of this effect have been recently made [s^, [tI, [tI, [tI, [tI, [t^, [ill ■ 
Another interesting direction is to develop a fast computation of non-linear power spectrum or correlation function 
for an arbitrary cosmological model. Recently, the statistical sampling method for precise power spectrum emulation 
has been proposed [zl, |80, HH . In this treatment, only a limited set of cosmological models can be used to predict 
power spectrum at the required accuracy over the prior parameter ranges. The analytic approaches combining this 
method may provide a fast and reliable way to estimate the two-point statistics, and the development of this method 
would be valuable. 
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APPENDIX A: STANDARD PERTURBATION THEORY UP TO THE TWO-LOOP ORDER 

In this Appendix, we briefly summarize the standard PT and derive a set of perturbative solutions. Based on these 
solutions, we obtain the analytic expressions for power spectrum up to the two-loop order. 

As we mentioned in Sec. lIII Al standard PT is the straightforward expansion of the quantity ^a, and the perturbative 
solutions are obtained by order-by-order treatment of Eq. (|2.2p . In order to systematically derive the solutions, the 
Einstcin-de Sitter (EdS) approximation is often used in the literature. In the EdS approximation, the matrix Qab given 
by Eq. (|2.3p is replaced with the one in the EdS universe, i.e., ^l{ri) = 1 and / — dlnD/dlna — 1. This means that 
all the non-linear growth factors appearing in the higher-order solutions are expressed in terms of the linear growth 
factor D(t). Neglecting the contributions from the decaying mode, the resultant solution for $a is then expanded as 

$ ,(fc; r;) = e^^i^ (fc; r;) + e2''$(2) (fc; ^) + e3''$(3) (fc; r;) -f • • • , (Al) 

The solution for each order of perturbation is expressed as 

/ '^(2^)3(n-'i)" '^D(fc-fci fc„)^l"^(fci,--- ,fc„)^o(fci)---<5o(fc„), (A2) 

where 5^ is the initial density fleld which we assume Gaussian statistic. The function J-"" is the symmetrized kernel 

(n) 

of the n-th order solutions. The explicit expressions for the kernel Fa is obtained from the recursion relation, which 
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can be derived by substituting the expansion ^M^ with (K2\i into Eq. ((2?2)) (e.g., |3a. I47ll8a. l8a| ) : 
FW(fei) = (l, 1), 

n-l 

• • . ,fc„) - ai;) ^ 7,,rf(qi,q2)Fi™)(/ci, • • ■ , fc„)^^("-™)(fe„_™+i, ••■,&„) (A3) 



m— 1 

(n) 



with = '^i + ■ ■ ■ + fcm and ^2 = '''m+i + • • • + fen- Here, the matrix a^^^ is given by 



(n) _ 1 ( 2n+l 2 1 . 

" (2n + 3)(n-l) I 3 2n ' • ^ ^ 



(n) 

Note that the kernel Fa given above is not yet symmetric under the permutations of arguments, fci, • • ■ , A;„, and it 
should be symmetrized: 

H''^-^^ E Fi^Hki,--- (A5) 

permutations 

Using the perturbative solutions, the power spectrum defined by (|3.ip is expanded as 

P,,{k; r^) = e^" P^l'\k) + e'^ { P^^\k) + P^^\k) } + e^" { P^f (A:) + P^f (fc) + P^f (fc) } + ■... (A6) 

Here, the quantity pC™") implies the ensemble average obtained from the m-th and n-th order perturbative solutions. 
In the above expression, the first term at the right-hand side is the linear power spectrum, while the second and third 
terms proportional to the growth factors e^"^ and e^"^ are respectively the so-called one-loop and two-loop corrections. 
The explicit expressions for these corrections become (e.g., pGlls^) 

= UaUbP^ik) (A7) 

Pff (fc) = 2 j -0^:Fi'Hq,k-q)Ti'\q,k-q)Po{q)Poi\k-q\) (A8) 

Pil'\k) ^ 3Po(fc) I -0^ {.Fp)(fe,q,-q)+^f (fc,g,-g)} Po(g) (A9) 

P^^\k) = 9Poik) I ^^^(3)(fc,p,-p)^f (fc,q,-q)Po(p)Po(g) 

+6 J ^^^H^\p,q,k~p^q)Ti'\p,q,k-p-q)Po{p)Po{q)Po{\k-p-q\) (AlO) 

Pif(fc) = 12 lf0^[T(')ip,k-p)Ti'\p,q,-q,k~p)+Ti'Hp,q,-q,k^p)^^^^^ 

X Po{p) Po{q) Po{\k - p\) (All) 
Pif (fc) = 15Po(A:) I ^^^{^(5)(p,q,fc,_p,-g) (p,q,fc,-p,-q)} Po(p)Po(<z), (A12) 

where Pq is the initial power spectrum of the density field Sq defined by Eq. ()2.7|) . and we set Ua = (1,1). 

Note that the expression for one-loop power spectra can be further reduced to the one-dimensional and two- 
dimensional integral for P^^^' and P^^^\ respectively (e.g., [H, [H, [s^ [s^ ) . In the results presented in Sec. IIVBTI 
we used the method of Gaussian quadratures for numerical integration of one-loop power spectra. On the other hand, 
for the two-loop power spectra, the integration cannot be simplified except for the first term in P^^^\ and we need 
to directly evaluate the six-dimensional integration. We adopted the Monte-Carlo integration to the two-loop power 
spectra. The integration kernels for each term are generated numerically using the recursion relation (jA4[) and the 
condition (IA5D. 



APPENDIX B: COMPARISON TO OTHER WORKS 



In this Appendix, we collect several recent works that attempt to improve the prediction of power spectrum and/or 
two-point correlation function, and discuss their qualitative differences. A quantitative aspect of various analytic 




FIG. 13: Diagrammatic representation for the perturbative treatment of the power spectrum proposed by Crocce & Scoccimarro 
(2008), based on the renormalized PT. This is compared with Fig. [3] 
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FIG. 14: Diagrammatic representation of the non-perturbative treatment proposed by Pietroni (2008), which can be compared 
with Fig. [1 



methods has been recently investigated in Ref. [4^. Here, we specifically comment on the approaches proposed by 
Refs. [13, m, [H, , which are very close to our treatment. 

Crocce &: Scoccimarro (2008) : First let us mention the work by Ref. 34]. Although the treatment presented 
in the paper are often quoted as 'RPT', strictly speaking, this is just the approximate treatment, which differs 
from the renormalized PT [32]. As we mentioned in Sec. IIII Al renormalized PT is the exact non-perturbative 
formulation without any approximations, and the power spectrum given by Eq. (j3.4p is expressed as the infinite 
series of irreducible loop diagrams constructed from the non-linear propagator, full vertex, and non-linear power 
spectrum. To make the analysis tractable, they adopted the following approximations: (i) the renormalized 
vertex is well-described by the (linear) vertex function; (ii) the non-linear power spectra that enter into the 

calculation of P^*^^ are all replaced with the linear-order ones. In our language, this corresponds to the first- 
order Born approximation. Then, using the approximate solution for propagator in Ref. [33|, they explicitly 
calculated the power spectrum including the corrections up to the two-loop order. The diagrams that they 
actually computed are shown in Fig. 1131 

Compared to our analytical treatment with Born approximation, there are two main differences. One is the 
higher-order corrections that appear in the diagrams (see Fig. [3]) . Another important difference is the asymptotic 
behaviors in the non-linear propagator. At ^ oo, the propagator used in their paper behaves like Gab — > 
(7Qf,exp[— a;^/2], which contrasts with Gab — > gab Ji{'2'X) / x in our closure approximation, where gab is the linear 
propagator and x is defined by a; = ka^{e'^ — e^). These distinctive features come from the partial resummation 
of a different class of higher-order terms when constructing the approximate solution of non-linear propagator 
(see Ref. [H, in details). Despite these remarkable differences, it has been shown in Ref. [i^ that the 
leading-order calculations neglecting the higher-order terms (two-loop diagram or second-order Born correction) 
can produce the same results which is indistinguishable from each other. This is true at least on large scales, 
where the agreement between N-body simulations and improved PT predictions is better than a few percent. 

Pietroni (2008) : Next consider the method proposed by Ref. [s^, called time-RG method. This method is based on 
the moment-based approach, and we first write down the moment equations. In general, this produces an infinite 
hierarchy of equations, however, Ref. [3^ assumes a vanishing trispectrum in order to truncate the hierarchy. 
As a result, a closed set of equations for power spectrum Pab is obtained, which coulples with the evolution of 
bispectrum Babe in some non-perturbative ways. Diagrammatic representation of this closed equations is shown 
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in Fig. [Ml which can be compared with Fig. [2] in our treatment. Note that in the subject of statistical theory 
of turbulence, this truncation procedure is referred to as quasi-normal approximation (e.g. .(4^. [ssl. [sojl. and it 
is known to have several drawbacks; positivity of the energy spectrum is not ensured, and it fails to recover the 
Kolmogolov spectrum in the inertial range of turbulence. 

Nevertheless, the advantage of this treatment, similar to our closure approximation, is that the power spectrum 
can be computed numerically by solving the evolution equations. This forward treatment seems quite efficient 
to bring out the non-perturbative effects incorporated into the formalism, and it has a wide applicability to 
include various physical effects. Recently, the formalism has been extended to deal with the effect of massive 
neutrinos [90l |. 

Valageas (2007) : The method proposed by Ref. [12] is based on the path-integral formalism. Starting from the 
action for the cosmological fluid equation (|2.2p . which describes the statistical properties of the vector field 
the large-N expansions as a technique of quantum field theory have been applied to derive the governing 
equations for power spectrum and propagator. In Ref. [i^, two kinds of expansions have been presented, leading 
to the two different non-perturbative schemes, i.e., steepest descent method and 2PI effective action method. 
Although both methods consistently reproduce the standard PT at the one-loop level, the latter includes the 
non-purturbative contributions which are not properly taken into account by the former method. Thus, the 
2PI effective action method is expected to provide a better result. It is interesting to note that despite the 
field-theoretical derivation, the resultant governing equations for the 2PI effective action method turn out to 
be mathematically equivalent to those obtained from the closure approximation [43| . Hence, the diagrammatic 
representation of this formalism is exactly the same as shown in Fig. [2] 

Matsubara (2008a) : Finally, we briefly mention the treatment proposed by Ref. [IHl- This is the Lagrangian- 
based approach, and we begin by writing down the exact expressions for matter power spectrum in terms of the 
displacement vectors. The resultant expression is in the exponential form, and the purterbative expansions are 
then applied for the explicit calculation of the ensemble average. While a naive expansion of the displacement 
vectors, together with the solutions of Lagrangian PT, merely reproduces the (standard) Eulerian PT results, 
Ref. |35j has applied a partial expansion, and some of the terms have been kept in the exponential form. This 
can be interpreted as the partial resummation of a class of the infinite diagrams. The resultant expression for 
power spectrum is quite similar to the one-loop result of standard PT, but slightly differs from it in the sense 
that there appears the exponential prefactor. As a consequence, the prediction reasonably recovers the damping 
behavior of the BAOs seen in the N-body simulations, and it also explains the smearing effect on the baryon 
acoustic peak in the two-point correlation function. 

One noticiable point of this method is that it is rather straightforward to generalize the calculations in real 
space to those in redshift space, since the displacement vectors in redshift space can be simply given by a linear 
mapping from those in real space. Further, the computational cost is less expensive compared to the other 
analytic methods. Although the validity range of this method is restricted to a narrow range of the low-fc 
modes, it would be very powerful for a fast compuation of the two-point correlation function. 

APPENDIX C: CONVERGENCE OF DIFFERENT COMPUTATIONAL METHODS FOR TWO-POINT 

CORRELATION FUNCTION 

In this paper, the grid-based calculation with FFT has been used for computing the two-point correlation functions 
from N-body data. Here, we compare it with other computational methods and check their convergences. 

Fig. 1151 shows the two-point correlation functions measured at z = 0.5 from the single realization of wmap5 simula- 
tions. Upper-left panel shows the results from the direct pair-counting. For each particle, we randomly select pairs, 
which are accumulated for each bin of separations, allowing for oversampling. The estimated values of two-point 
correlation function are then plotted for different number of samples A^sampio: A^sampio = 640k, 2, 560k, 10, 240k and 
40, 960k. The resultant total number of pairs, A^pair, indicated in the panel is given by Apair = Agampic x A'particic, with 
A^particio = 512'^ being the total number of particles. Note that the actual number of pairs that enters into the plotted 
range is less than A'pair. On the other hand, jrpper-right panel shows the results from the grid-based pair-counting 
introduced by Barriga & Gaztaiiaga (2002) [91| (see also Ref. [Z^)- this method, we first construct the density 
field on a grid of A^grid cells, and then estimate the correlation function through the pair count on grids: 
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FIG. 15: Comparison between different computational methods for two-point correlation function. 



Compared to the direct pair-counting, this method is computationally efficient when we store the list of neighbor 
particles which contribute to a given bin of separation. We plot the results adopting the two different number of cells, 
-^grid — 128"^ and 256'^. In lower-left panel, the grid-based calculation with FFT (see Eq. (|4.2p ) is used to compute 
the two-point correlation function, with different numbers of cells, iVgiid = 128"^, 256"^, 512"^ and 1, 024"^. Note that we 
adopt A^grid = 1,024^ in the analysis presented in Sec. IIVI Finally, in lower- right panel, the results for three different 
methods with the largest number of pairs or grids are collected and compared with each other. 

To check the convergence, we further evaluate the residuals from the mean values, = ^ — ^, and plot the results 
in each panel of Fig. 1151 Here, the mean values ^ are estimated from the ensemble average over the three different 
results using the largest number of pairs or grids. As increasing the numbers A^pair or A'grid, the results for three 
different methods all approach the mean values ^, and a few percent-level agreement is achieved over the range of 
our interest (except for the vicinity of zero-crossing point, ^ w 0). It is interesting to note that residuals obtained 
from the grid-based pair-count and FFT methods almost coincide with each other and the differences are hard to 
distinguish, indicating that both methods are equivalent even in the practical situation. These experiments suggest 
that the grid-based calculation with FFT is a reliable estimation method comparable to the other methods. It should 
be emphasized that the method using FFT is much more efficient than other pair-count methods. For example, using 
8 cores of 3GHz processors, the direct pair-counting with iVpair = 10, 240k x A'particie takes about two weeks to get the 
results shown in Fig. [151 The grid-based pair-counting is computationally less expensive than the direct pair-counting, 
but it still needs time-consuming calculations, especially for a large number of grids. By contrast, the method using 
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FFT only requires few minutes even with Ng^id = 1,024^. This can be achieved by a single- node calculation. 
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